Optical properties of bulk semiconductors and graphene/boron-nitride: The 
Bethe-Salpeter equation with derivative discontinuity-corrected DFT energies 
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We present an efficient implementation ol the Bethe-Salpeter equation (BSE) for optical prop- 
erties of materials in the projector augmented wave method GPAW. Single-particle energies and 
wave functions are obtained from the GhLBSC functional which explicitly includes the derivative 
discontinuity, is computationally inexpensive, and yields excellent fundamental gaps. Electron-hole 
interactions are included through the BSE using the statically screened interaction evaluated in 
the random phase approximation. For a representative set of semiconductors and insulators we find 
excellent agreement with experiments for the dielectric functions, onset of absorption, and lowest ex- 
citonic features. For the two-dimensional systems of graphene and hexagonal boron-nitride (h-BN) 
we find good agreement with previous many-body calculations. For the graphene/h-BN interface, 
we find that the fundamental and optical gaps of the h-BN layer are reduced by 2.0 eV and 0.7 eV, 
respectively, compared to freestanding h-BN. This reduction is due to image charge screening which 
shows up in the GLLBSC calculation as a reduction (vanishing) of the derivative discontinuity. 

PACS numbers: 71.15.-m, 78.20.-e, 71.35.Cc 
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I. INTRODUCTION 

Optical spectroscopies such as photo absorption, lumi- 
nescence, and reflectance measurements are widely used 
for materials characterization. In this context, first- 
principles calculations play an increasingly important 
role for the interpretation and guidance of experimental 
investigations. However, theoretical spectroscopic meth- 
ods are not only useful for characterization purposes. 
Indeed, with the recent focus on solar energy conver- 
sion, plasmonics, and optoelectronics - all applications 
which involve the interaction of light with matter - first- 
principles methods for calculating the optical properties 
of complex materials are becoming the essential tool al- 
lowing for reliable computational design of new materials 
within these areas. 

The two most commonly used ab-initio methods for 
optical properties are time-dependent density functional 
theory (TDDFT)A and many-body perturbation theory 
(MBPT) 2 . For smaller molecules and clusters 3 -, TDDFT 
with the adiabatic local density approximation (ALDA) 
provides a reasonably good compromise between accu- 
racy and computational cost. However, the ALDA fails 
to describe several important effects including the for- 
mation of excitons in extended systems^, charge-transfer 
excitations in donor-acceptor molecular complexes^, as 
well as the screening of optical transitions by nearby 
metal surfaces^. Apart from these qualitative failures, 
the ALDA is also found to underestimate the optical 
transition energies and overestimate static dielectric con- 
stants of bulk insulators and semiconductors. This prob- 
lem is, at least to some extent, related to the well known 



tendency of the LDA and related semi local exchange- 
correlation (xc) functionals, to underestimate the funda- 
mental energy gaps in such systems. 

All of the above mentioned problems of the TDDFT- 
ALDA approach are overcomed by the MBPT. In the 
standard scheme, the quasiparticle band structures are 
obtained using the GW approximation^ while optical ex- 
citation energies are obtained by solving a Bethe-Salpeter 
equation (BSE)£ with a statically screened electron-hole 
interaction. The GW-BSE approaches has been sue- 
cesfully applied to a number of different systems rang- 
ing from bulk semiconductors 9 , insulators and their 
surfaces^, two-dimensional systems such as graphene 12 
and boron nitride layers^ 3 -, metal-molecule interfaces*, 
isolated molecules^— and liquid water—. Nevertheless, 
applications of the approach to larger systems are limited 
by the extremely demanding computational requirements 
of both the GW and BSE calculations. 

Several schemes have been proposed to reduce the 
computational cost of GW-BSE calculations. These 
include circumventing the GW step by applying sim- 
pler band structures e.g. derived from the COHSEX 
approximation^ or simply scissors operator-corrected 
LDA band structures^ 9 -, or the use of model dielectric 
functions to describe the screening 2 ^. Another route of 
research is directed towards the development of more ac- 
curate TDDFT xc-kernels without sacrificing the compu- 
tational simplicity associated with this approach^—. 

Recently, Kuisma et al. have introduced the GLLBSC 
xc-potential2^ which is based on an earlier functional de- 
veloped by Gritsenko et al£^. This potential explicitly 
includes the derivative discontinuity of the xc-potential 
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at integer particle numbers which is important to obtain 
physically meaningful band gaps from DFT. The deriva- 
tive discontinuity, A xc , is calculated directly from the 
Kohn-Sham eigenvalues and eigenstates. The fundamen- 
tal band gap is then obtained as the sum of the Kohn- 
Sham single-particle gap and the derivative discontinu- 
ity. The GLLBSC method has been shown to produce 
fundamental band gaps as well as band dispersions for a 
range of semiconductors in very good agreement with ex- 
periments and more sophisticated theoretical approaches 
while the computational cost is comparable to that of 



In this paper we combine the TDDFT and BSE meth- 
ods for treating the electron-hole interaction with the 
GLLBSC method for the wave function and band struc- 
tures. Considering both bulk and low dimensional sys- 
tems we find that the accuracy of the GLLBSC-BSE ap- 
proach is comparable to the GW-BSE approach. All 
the methods are implemented in the GPAW code22r— , 
an electronic structure package based on the projector 
augmented wave methodolog y 31 ! 32 . For the bulk sys- 
tems Si, C, InP, MgO, GaAs and LiF, we find that 
the fundamental gaps and static dielectric constants cal- 
culated with GLLBSC compare well with experimental 
data. Importantly, the static dielectric constant should 
be evaluated without the derivative discontinuity when 
using an xc-kernel that does not account for e-h interac- 
tion such as the ALDA or the random phase approxima- 
tion (RPA). The experimental optical absorption spec- 
tra of all compounds are also very well reproduced by 
the GLLBSC-BSE approach including the absorption on- 
set and excitonic peaks. Finally, the method is used to 
compute the band structure and optical absorption spec- 
tra of graphene, hexagonal boron- nitride (h-BN), and a 
graphene/h-BN interface. For the isolated sheets we find 
good agreement with previous GW-BSE calculations. 
For the interface we find that both the quasiparticle- and 
optical gap of the h-BN sheet are reduced by 2.0 and 0.7 
eV, respectively. The physical origin of this effect is due 
to image charge screening by the graphene layer. In the 
GLLBSC, the reduction shows up as a vanishing of the 
derivative discontinuity. 



The rest of the paper is organized as follows. Section 
II introduces the theoretical framework for calculating 
optical properties of solids with GPAW using the TDDFT 
and BSE approaches, followed by a brief review of the 
GLLBSC method. Details of the implementation are pre- 
sented in Sec. III. Section IV presents benchmark results 
for the band gaps, dielectric constants and optical ab- 
sorption spectra of a number of bulk semiconductors and 
insulators. In Sec. V we present the band structures and 
optical spectra of graphene, h-BN, and graphene/h-BN 
interface. Finally, a summary is given in Sec. V. 



II. METHOD 
A. Macroscopic dielectric function 

Most of the optical properties of a solid can be obtained 
from the macroscopic dielectric function, 



(i) 



G=0,G'=0 



Here, £gg' (q> w ) is the (microscopic) dielectric matrix in 
reciprocal G space. The off-diagonal elements of the e 
matrix account for local field effects arising due to the 
periodic crystal potential. The macroscopic average is 
achieved through the inversion of the e matrix. 

In this work we consider only the longitudinal compo- 
nent of the dielectric function. For applications to optical 
properties this is in fact not a restriction because in the 
relevant long wave length limit the electrons do not feel 
the difference between longitudinal and transversely po- 
larized fields, and consequently the two types of response 
functions coincide. Still, for anisotropic systems e(ui) de- 
pends on the direction in which the limit q — ¥ is taken. 
However, to keep the notation simple we shall omit ref- 
erence to this direction in what follows. 



B. Linear response function from TDDFT 

The microscopic dielectric matrix is related to the lin- 
ear density response function, \, via 
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'GG' 



|q+G||q + G 



7TXGG'(q,w). (2) 



Within TDDFT the response function is related to the 
response function of the non-interacting Kohn-Sham elec- 
trons, x° an d the exchange-correlation interaction kernel 
K xc via a Dyson-like equation, 

XGc(q,w) = x GG '(q,w) 

+ E Gl G 2 X% Gl (q, u)K Gi g 2 (q, w)xg 2 c (q, £*>)• (3) 



The KS response function is given by 
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GG'(q^) = 7^ X! ^ nk ~~ /«'k+q) 



k,)m' 



lnk,n'k+q(G)n,* k ?l , k+q (G') 
W + £„k - En'k+q + if} 



(4) 



where e n \^ is a KS eigenvalue, and /„k is the occupation 
factor. The quantity 



«nk,n<k+ q (G) = (V>„k|e- l(q+G) - r |Wk +q ) 



(5) 



is referred to as the charge density matrix^.. In the long 
wavelength limit, i.e. for q — > 0, and for n 7^ n' , applica- 
tion of the k-p perturbation theory 35 yields the important 
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identity 



lim q ^O«nk,r!'k+q(0) 



-^q • (^KklviV'n'k) 

^n'k ^nk 



(6) 



Alternatively, this form follows directly if we consider the 
density induced by a longitudinal vector potential rather 
than a scalar potential. A detailed description of the 
evaluation of the charge density matrix and the ALDA 
xc-kcrncl within the PAW formalism can be found in Ref. 

M 



C. The Bethe-Salpeter Equation 

Several of the shortcomings of the ALDA in describing 
optical spectra are overcomed by explicitly accounting 
for electron self-energy effects and electron-hole interac- 
tions using many-body perturbation theory. In the stan- 
dard GW-BSE approach, the single-particle energies are 
evaluated using a self-energy in the GW approximation 
while the optical excitation energies are obtained by di- 
agonalizing an effective two-particle Hamiltonian. In the 
present work we avoid calculating the GW self-energy by 
using single-particle energies obtained from the efficient 
GLLBSC functional. 

Following the standard approach, the excitation ener- 
gies corresponding to an external potential with momen- 
tum q can be found by solving an eigenvalue problem of 
the form 



]>>(q)ss^,(q)=£: A (q)4(q) 



(7) 



where 'Hss'(q) is the Bethe-Salpeter effective two- 
particle Hamiltonian evaluated in a basis of electron- 
hole states, Vs(rfc,r e ) = -0nk( r /O* ? /'mk+q(r e ). The BSE 
Hamiltonian reads 



UsS'iq) = (e^k+q - £ nk ) S SS' ~ (/mk+q - fnk)K SS ,(q) 



The kernel consists of an e-h exchange interaction (V) 
and a direct screened e-h attraction (W), 



^ss'(q) =V SS '(q)- ^ss'Cq) 



(9) 



The factor 2 accounts for spin. In appendix |A1 we give a 
derivation of the BSE eigenvalue equation and its relation 
to the dielectric function. 

The effective two particle Hamiltonian is most conve- 
niently evaluated in a plane wave basis. In this represen- 
tation the e-h exchange term reads 



VW(q) = 



4tt 



Tl ?lk,?7lk+q 



(G)n r , 



-q(G) 



G 



|q + GP 



(10) 



If we exclude the G = component in the sum we ob- 
tain the short range exchange kernel V. The difference 
between V and V becomes important when the response 



function is written in terms of the eigenstates and ener- 
gies of the BSE Hamiltionian, see below. To obtain the 
optical limit Vss' (q — > 0) we use the expression Eq. ([6]) 
to cancel the 1/q 2 Coulomb divergence appearing in the 
G = term. In the evaluation of the remaining terms 
we use a small finite value for q (a value of 0.0001 A -1 
has been used in this work). 

The plane wave expression for the e-h direct Coulomb 
term reads 



4tt 



W SS ,(q) = — ]T < k ,„, k ,(G)W GG >(k' - k) 



GG' 



X n mk+q,m'k'+q(G'), 



where 

W GG <(k'-k) 



GG' 



(k'-k,w = 0) 



Ik'-k + Gllk'-k + G'l 



(11) 



(12) 



Here we encounter a divergence of W GG ' when either G 
or G' is zero and k = k'. Such a divergence due to the 
singularity of the Coulomb kernel at q = is also present 
in calculating exact exchange^ and GW self energies^. 
When n ^ n' and m ^ m' we can use the expression Eq. 
© to cancel the divergence; while for n = n' or m = mf, 
the singularity in the Coulomb kernel is integrated out 
analytically, following Ref. around a sphere centered 
at q = 0. We have also adopted another scheme using 
an auxiliary periodic function with the same singular- 
ity as the exact function but which can be evaluated 
analytically^. These two schemes give essentially the 
same results. 

The eigenstates and eigenvalues of the BSE Hamilto- 
nian provide a spectral representation of the four-point 
density response function (see Appendix |A")) . 



4P , , A _ V ^(q)[^(q)mv 

xss'iq^J - 2^ 



XX' 



E x (q) + ir\ 



(13) 



(8) where N\\> is the overlap matrix defined as 



jvju, = 5>5(q)]*4&'(q). (i4) 

s 

Using the plane wave representation ([5]) of the electron- 
hole basis states we obtain the following expression for 
the response function in reciprocal space 

XGG'(q,w) - ^J> 4 5 P 5<(q^)"s(G)n*,(G') (15) 



SS' 



From this expression the inverse dielectric constant and 
macroscopic dielectric constant follows from Eq. ([2]) and 
(HJ), respectively. 

We note that upon excluding the 1/q 2 term in the e-h 
exchange term, i.e. replacing V by V in the kernel ([9]), 
the eigenstates and eigenvalues of the BSE Hamiltonian 
provides a spectral representation of the irreducible re- 
sponse function— rather than the full response function. 
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In this case the effect of V is to account for local field 
effects. Consequently the macroscopic dielectric function 
can be written 



4-7T 

e(w) = 1 - p jgXooiq -> 0, w) 



47T v - 



n s (0)n^(0)/ s 



AA' 



(16) 



In the above expression the optical limit q — >• is taken 
in the following way. First, the BSE Hamiltonian is con- 
structed using an e-h basis of vertical excitations (q = 0) 
but using a finite small q for the Coulomb interaction 
l/|q + G| in V (or V). The same finite q is then used 
when evaluating the dielectric function from the spectral 
representation of the (irreducible) response function. 



D. Quasiparticle energies from GLLBSC 

The derivative discontinuity A xc is defined as the dif- 
ference between the fundamental gap E s and the Kohn- 
Sham (KS) single-particle gap as follows 

E g = I- A = E{n N ^}-2E[n N } + E[n N+1 ] = E* s +A xc , 

(17) 

where E[tim\ is the total energy of the TV-electron system 
and the fundamental band gap E s is defined as the dif- 
ference betweeen the ionization energy I and the electron 
affinity A. 

Within the GLLBSC method, the derivative disconti- 
nuity A xc is obtained through 



(V N+1 \A(r)\V N+1 ) 



(18) 



where 



occ 



A( r ) = )■ y Kx [Vclumo - f-i ~ \Ahomo - £ij 



l^(r)| s 
n(r) 



(19) 

£i, ipi{r) and n(r) are eigenvalues, eigenstates and elec- 
tron density, respectively, obtained from solving the KS 
equation with the following GLLBSC potential 



«GLLBSc(r)=2 e P c BEsol (r) 



(20) 



n(r) 



c,rcsp 



Here, K x ss 0.382 is a coefficient fitted from electron gas 
calculations to reproduce the exchange potential for uni- 
form electron density and e r is a reference energy taken 
from the highest occupied eigenvalue. The GLLBSC 
method is an orbital dependent simplification of the KLI 
approximation to the exact-exchange optimized effective- 
potential method following the guidelines of GLLB 2 ^ for 
the exchange potential. For the details of the formulation 
we refer the reader to Ref. [24. 



III. IMPLEMENTATION 

The TDDFT and BSE codes are implemented in 
GPAw22r— , a real-space electronic structure code using 
the projector augmented wave methodolog y 31 ' 32 . In this 
section, we focus on the construction of the screened 
Coulomb interaction kernel W, which is the most chal- 
lenging and time consuming part in the BSE formalism. 
For the details of the implementation on the GLLBSC 
potential and the linear density response function in the 
PAW formalism, we refer to Ref. |24| and l30l. respectively. 



A. Screened Coulomb interaction W 

The electron-hole correlation kernel Eq. (ITU contains 
the dynamically screened Coulomb interaction in a plane 
wave representation, 



w G g' (q, w ) = 



te G g (q, u) 
q+G||q + G' 



(21) 



In Eq. (jlll) the q vector represents the difference between 
two k-points in the first Brillouin zone. Thus, the q-point 
mesh has the same form as the k-point mesh. In addition, 
the q-point mesh always includes the T point, while the 
k-point mesh does not necessarily. The use of k-point 
symmetry for obtaining the wave-functions at k-points 
outside the irreducible Brillouin zone has been described 
in a previous paper—. In the following we describe how 
symmetry considerations can be used to reduce the q- 
point sum. 

We start by examining the q-point symmetry in the 
charge density matrix defined in Eq. ([5]). Consider a q 
satisfying 



q = T qmz + G o 



(22) 



where q IBZ is an irreducible q point, T is a crystal sym- 
metry transformation, and Go is a reciprocal lattice vec- 
tor that translates the Tq IBZ vector back into the Bril- 
louin zone if needed. The charge density matrix in Eq. 
([5]) then becomes 

^nk,n'k+q(G) 

= (Vnk|e- i ^- +G ° +G > r |Vn'k +q ) 

= (^T- 1 k|e- j[q - +T " (Go+G)I - r |^T- 1(k+q )> 

= ^nT- 1 k,n'T- 1 (k+q)(T _1 (G + G)) (23) 

Since the calculation of Xgg' (l> w ) involves the sum- 
mation of the charge density matrix over all the BZ k- 
points, the above equation leads directly to the following 
relation (as long as T _1 k belongs to the k-point mesh): 

XGG'(q> w ) = XT-i(G+Go),T-i(G'+G )(qiBZ>w). (24) 

The above relation also applies to W / GG'(q> a; )- 
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Besides crystal symmetry, time reversal symmetry is 
also used for systems that have no inversion symmetry. 
If the transformation of a given q to IBZ requires both 
crystal symmetry and time reversal symmetry via 



(25) 



q = -^qiBz + G , 
the W matrix should satisfy 

WbG'fa.w) = W-T-HG+Go),-T-i(G'+G )( < llBZ^)' 

(26) 

Finally, it has to be emphasized that for a hnte k-point 
mesh used in a numerical calculation, the crystal symme- 
try transformation T should apply to both q-points and 
k-points. This results in reduced crystal symmetry oper- 
ations if the r centered q-point mesh does not coincide 
with the k-point mesh. 



IV. SOLIDS 

In this section the optical properties of a representative 
set of six bulk semiconductors and insulators are studied 
using both ALDA and the BSE. We start by presenting 
the fundamental gaps obtained with LDA and GLLBSC. 
The accuracy of the GLLBSC gaps is similar to Go Wo 
calculations from the litterature with an average abso- 
lute deviation of 0.3 eV from experiments. An impor- 
tant ingredient in the BSE calculation of optical spec- 
tra is the static dielectric constant which determines the 
strength of the screened electron-hole interaction, W. We 
find that the best agreement with experiment is obtained 
when the response function is evaluated from the LDA 
or GLLBSC Kohn-Sham (i.e. without adding the deriva- 
tive discontinuity) energies, and we explain this from the 
fact that the electron-hole interaction is not explicitly 
accounted for by the random phase approximation used 
to obtain e. Finally, the absorption spectra using both 
ALDA and BSE are presented. Very good agreement 
with the experimental spectra is found for the GLLBSC- 
BSE combination both for the absorption onset and the 
excitonic features. 



A. Fundamental gaps 

Table I shows the calculated band gaps for Si, C, InP, 
MgO, GaAs and LiF. We have used the experimental 
lattice constants for all systems: Si (5.431 A), C (3.567 
A), InP (5.869 A), MgO (4.212 A), GaAs (5.650 A) and 
LiF (4.024 A). The Kohn-Sham energies and wave func- 
tions were obtained with GPAW using uniform grids with 
spacing 0.2 A and a Fermi temperature of 0.001 eV. The 
Brillouin zone was sampled using a Monkhorst-Pack grid 
of 24 x 24 x 24 which was found sufficient to converge the 
band gaps to within 0.02 eV. 

Compared to LDA band gaps (first column), GLLBSC 
even without the discontinuity (second column) improves 



TABLE I: Band gaps (units in eV) calculated using GLLBSC 
without (wo.) and with (w.) the derivative discontinuity A xc 
added to the Kohn-Sham gap. These values are compared 
with LDA, Go Wo and experimental data. Underlined values 
correspond to zero-temperature values. The mean absolute 
errors (MAE) with respect to experiments are summarized in 
the last row. 





LDA 


GLLBSC 
(wo.) 


GLLBSC 
(w.) 


G Wo 


Expt. 


Si 


0.51 


0.74 


1.09 


1.12 a 


1.17 b 


c 


4.16 


4.22 


5.52 


5.50" 


5.48 c 


InP 


0.61 


1.15 


1.63 


1.32 d 


L42 6 


MgO 


4.63 


6.10 


8.32 


7.25" 


7.83 e 


GaAs 


0.57 


0.93 


1.23 


1.30" 


1.52 b 


LiF 


8.87 


10.97 


14.94 


13.27" 


14.20 / 


MAE 


2.04 


1.25 


0.31 


0.32 





"Reference 3i| 
'Reference 4(1 
c Reference 4l] 
d Reference 43 
e Reference 431 
^Reference 44l 



T=0K 



TABLE II: The static macroscopic dielectric constant e ob- 
tained using TDDFT on top of LDA as well as GLLBSC elec- 
tronic structure without (wo.) and with (w.) discontinuity 
A xc applied. The two rows for each semiconductor correspond 
to TDDFT calculations with RPA and the ALDA kernel, re- 
spectively. 





LDA 


GLLBSC 
(wo.) 


GLLBSC 
(w.) 


Expt. 


Si (RPA) 


12.53 


11.00 


10.25 


11.90 9 


(ALDA) 


13.16 


11.54 


10.73 




C 


5.56 


5.48 


5.04 


5.70 9 




5.82 


5.74 


5.25 




InP 


11.48 


8.92 


8.06 


12. 5 9 




11.99 


9.33 


8.41 




MgO 


3.06 


2.52 


2.31 


2.95 ft 




3.20 


2.63 


2.39 




GaAs 


13.52 


11.12 


10.28 


11. 10 9 




14.17 


11.68 


10.78 





»Reference[47l T=300K. 

^Reference 4£| optical dielectric constant. 



the band gaps. The reason is that the GLLBSC po- 
tential Eq. (f2"0")) can reproduce the asymptotic 1/r be- 
havior of the Coulomb potential^ and thus the Kohn- 
Sham eigenvalues are improved over LDA. By adding 
the discontinuity (third column), the band gaps agree 
reasonably well with experimental data (last column). 
The mean absolute error (MAE) with respect to the ex- 
perimental data is 0.31 eV in agreement with a previ- 
ous study using GLLBSC for oxides in the perovskite 
structure^. The sign of the deviations from experiment 
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FIG. 1: Optical absorption spectra calculated using LDA-ALDA (dash-dotted line), GLLBSC-ALDA (dashed line) as well as 
GLLBSC-BSE (solid line). The derivative discontinuity, A xc , is included in the GLLBSC calculations. The calculated spectra 
are compared with experimental data (dots, Ref. |49| ). 



seem to vary randomly. This is in contrast to the Go Wo 
results (fourth column)^, which systematically underes- 
timates the band gaps with the largest error being almost 
1 eV. We note that (quasi-) selfconsistent GW calcula- 
tions have been shown to improve the ionization poten- 
tials of molecules^ and band gaps of solids^ by reducing 
the overscreening resulting from the LDA starting point. 
However, such calculations are even more computation- 
ally demanding than Go Wo, and are therefore not nor- 
mally used for the calculation of optical spectra. We will 
show in the following that GLLBSC represents a cheap 
alternative means to GW providing not only reasonable 
fundamental gaps, but also very good optical dielectric 
constants and absorption spectra. 

B. Dielectric constants 

Table II shows the calculated static macroscopic di- 
electric constants. In addition to the parameters pre- 
sented for obtaining the band gaps, 60 - 90 unoccupied 
bands, corresponding to around 140 eV above the Fermi 
level, were used in the calculation of the response func- 
tion Eq. @. Local field effects were included up to 
an energy cutoff of 150 - 250 eV, which varies accord- 
ing to the size of the unit cell and corresponds to 169 
G vectors. The static dielectric constants obtained us- 
ing LDA-RPA (first column), that is, RPA calculations 
based on LDA wave functions and energies, are generally 



higher than the experimental values (last column) due to 
the underestimated LDA band gaps. The overestimation 
is enhanced by inclusion of the ALDA kernel (the second 
row for each semiconductor) , in agreement with previous 
studies^. The GLLBSC without the discontinuity in- 
crease the band gaps relative to LDA and consequently 
reduces the dielectric function towards the experimental 
value. The inclusion of the discontinuity further opens up 
the gap and the corresponding dielectric constants (third 
column) systematically underestimate the experimental 
values. This underestimation is a result of the neglect 
of electron-hole interaction when the response function is 
evaluated at the RPA and (to some extent) ALDA levels. 
In order to reduce the error coming from this effect, the 
response function should be evaluated using "dressed" 
single-particle energies rather than the bare QP ener- 
gies. In the following, we use the GLLBSC(wo.)-RPA 
dielectric function for calculating W. 

C. Absorption spectra 

The absorption spectra calculated using TDDFT and 
the BSE are shown in Fig. 1. TDDFT calculations were 
performed using the ALDA kernel and the same param- 
eters as used for obtaining the dielectric constants (see 
previous section). For the BSE calculations we used an 
8x8x8 Monkhorst-Pack /c-point grid not containing the 
Gamma-point (for InP 10 x 10 x 10 fc-points were used). 
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FIG. 2: Band structure of a h-BN sheet calculated with 
GLLBSC (solid lines) and LDA (dotted lines). The top of 
the valence bands is set to zero. 

We have also checked the spectra with 12x12x12 fc-point 
sampling. The main peaks in the absorption spectra are 
well converged with the applied fc-point sampling, how- 
ever, a complete elimination of the small "wiggles" seen 
in the spectra would require significantly denser fc-point 
sampling. The screened interaction kernel, W"gg'(ci) 5 
was obtained using GLLBSC(wo.)-RPA, with 60 un- 
occupied bands and local field effects included by 169 
G- vectors. Three valence and three conduction bands 
were taken into account in constructing the BSE ma- 
trix. Again, this is sufficient to converge the major (ex- 
citonic) peaks and the low energy part of the absorption 
spectra. The Tamm-DancofF approximation^, consisting 
of the neglect of coupling between v-c and c-v transi- 
tions, was employed. The effect of temperature, which in 
general lowers the band gap and smears the absorption 
spectrum 5 -^, is not considered in the current work. As 
a result, the spectra presented here are broadened using 
smearing factors (in units of eV): Si (0.10), C (0.35), InP 
(0.20), MgO (0.25), GaAs (0.20) and LiF (0.12). 

As can be seen from the absorption spectra in Fig. [TJ 
LDA-ALDA (green dash-dotted lines) gives threshold op- 
tical transition energies that are 0.5-3 eV lower than 
experiments (black dots). This is a result of the too 
low LDA band gaps. The use of GLLBSC wave func- 
tions and energies including the derivative discontinuity, 
GLLBSC- ALDA (blue dashed lines) increases the absorp- 
tion threshold energies and improves the agreement with 
experiments. However, the shape of the spectra are qual- 
italitvely different. In particular, the spectra are too low 
at the on-set of the absorption and the excitonic features 
in Si, MgO, and LiF, are completely missed. This is be- 
cause ALDA does not properly account for electron-hole 
interactions. In contrast the spectra obtained from the 
BSE using the GLLBSC eigenvalues as QP energies (red 
lines) are in excellent agreement with experiments. A 
small exception is for GaAs where a small peak, abscent 



in the experimental spectrum, is seen at around 2 eV. 
A similar feature was seen in a previous calculaton em- 
ploying a non-local approximation to the xc kernel within 
TDDFT22, but does not appear in a previous GW-BSE 
calculation 51 . This indicates that the presence of the fea- 
ture is related to differences between the GLLBSC and 
GW band structure. We note that (small) deviations 
between the GLLBSC and GW band structures was re- 
cently proposed as the reason for (slight) inaccuracies in 
the GLLBSC- ALDA calculated surface plasmon energies 
ofAg(lll)21. 



V. GRAPHENE/BORON-NITRIDE 

In this section we study the bandstructure and optical 
absorption spectra of graphene, a single layer of hexago- 
nal boron-nitride (h-BN), and their interface graphene/h- 
BN. The lattice parameter of h-BN is very similar to that 
of graphene making it a promising candidate substrate 
material for graphene based devices 5 -^. In contrast to 
graphene, which is a semi-metal, h-BN has a wide band 
gap and exhibits strong excitonic effects. The optical 
properties of layered BN sheets as well as BN nanotubes 
have been studied extensively both experimentally 5 - 3 - and 
theoreticall y 13 ' 54 . Upon adsorption of graphene onto a 
h-BN sheet, a small bandgap of around 10-200 meV, 
depending on the configuration and interplane distance, 
emerges 55 . The ground state electronic properties, in- 
cluding the role of dispersive forces, and the band struc- 
ture have been studie d 55 ' 56 . Below we investigate the 
optical properties of the graphene/h-BN interface and 
assess the quality of the GLLBSC for such 2D structure. 

Before presenting the results for graphene/h-BN, first 
we examine a single h-BN sheet. For the lattice constant 
of h-BN we used 2.89A and 20A vacuum was included 
between the periodically repeated BN layers. Figure [2] 
shows the band structure calculated using LDA (dotted 
lines) and GLLBSC (solid lines). The LDA band gap (sit- 
uated at the K-point) is 4.61 eV which is 0.3 eV larger 
than reported in an earlier pseudopotential study—. The 
GLLBSC band gap is 7.99 eV, which includes the deriva- 
tive discontinuity of 2.12 eV, is close to the pseudopoten- 
tial G W band gap of 7.9 eV-51 

Figure [3] shows the absorption spectrum of a h-BN 
sheet obtained with three different methods. The LDA- 
ALDA spectrum shows a broad absorption peak with 
an onset at 4.5 eV in good agreement with literature 5 -^. 
The GLLBSC-ALDA spectrum is essentially identical 
to LDA-ALDA, but blue shifted by the difference in 
the band gap. For the BSE calculation, the Brillouin 
zone was sampled on a non Gamma-centered 32 x 32 
Monkhorst-Pack grid, and 70 unoccupied bands were in- 
cluded to obtain the screened interaction W . A two- 
dimensional Coulomb cutoff technique^ was used to 
avoid interactions between supercells. Since we are inter- 
ested in the low-energy part of the absorption spectrum 
and because the valence and conduction bands are well 
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FIG. 3: Optical absorption spectra of a h-BN sheet calcu- 
lated using LDA-ALDA (dash-dotted line), GLLBSC-ALDA 
(dashed line) and GLLBSC-BSE (solid line). 
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FIG. 4: Top (a) and side (b) view of a graphene/h-BN. (c) 
Band structure of graphene/h-BN calculated with GLLBSC 
(solid lines) and LDA (dotted lines). The top of the valence 
bands is set to zero. 



separated from the rest of the bands in the relevant part 
of the Brillouin zone (around the K-point), only the va- 
lence and conduction bands were included in the BSE ef- 
fective Hamiltonian. The absorption spectrum obtained 
with GLLBSC-BSE shows three excitonic peaks at 6.1, 
7.1 and 7.4 eV with decreasing amplitude. These exciton 
energies agree well with the value of 6.2 eV, 7.0 eV and 
7.4 eV obtained with the GW-BSE scheme^. 

For the graphene/h-BN interface, we studied the struc- 
ture where one C atom is ontop of a B atom and the other 
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FIG. 5: Upper panel: Optical absorption spectrum of 
graphene/h-BN calculated using LDA-ALDA (dash-dotted 
line), GLLBSC-ALDA (dashed line) and GLLBSC-BSE (solid 
line). Lower panel: The GLLBSC-BSE spectrum of the in- 
terface (repeated) together with the sum of the absorption 
spectra of an isolated graphene and BN layer, respectively. 



C atom is above the center of the BN ring, as shown in 
Fig. [4] (a) and (b). Both graphene and h-BN are kept 
planar at a distance 3.48A apart. Recent RPA calcu- 
lations found this structure and adsorption distance to 
be the most stable^. Fig. @] shows the band structure 
of graphene/h-BN. For the LDA band structure (dotted 
lines), a small band gap of 31 meV opens at the K point. 
This number is very close to the 53 meV found in an ear- 
lier study^. The h-BN gap, indicated by the arrow and 
is 4.60 eV in the LDA, which is essentially the same as 
found for the isolated h-BN sheet (4.61 eV). This is in 
contrast to the GLLBSC band structure which yields a 
band gap of the adsorbed h-BN of 6.01 eV which is 1.98 
eV lower than obtained for isolated h-BN. This sizable re- 
duction of the gap is not due to hybridization, but rather 
is a result of a reduction of the derivative discontinuity 
from 2.12 eV to essentially zero. We note in passing that 
the GLLBSC value of 6.01 eV is 0.3 eV larger than our 
Go Wo results for this system (to be published elsewhere). 

The reduction of the fundamental gap when BN is 
adsorbed on graphene is physically meaningful and can 
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be explained by the screening provided by the graphene 
layer (image charge effect) which reduces the energy cost 
of removing electrons/holes from the BN layer. For 
molecules on surfaces, this effect has been shown to be 
well described by the GW method, whereas both (semi- 
)local and hybrid functionals completely miss the effect 
predicting no change in the gap upon adsorption (apart 
from obvious hybridization effects) 59 ' 60 . Interestingly, 
within the GLLBSC the gap reduction is a result of the 
vanishing, or strong reduction, of the derivative discon- 
tinuity. However, this also has the unphysical conse- 
quence that the reduction is present independent of the 
graphene-BN distance. This follows from the obsrvation 
that the derivative discontinuity in Eq. (fT"9|) becomes 
zero for a metallic system. 

The absorption spectrum of graphene/BN calculated 
with the three different schemes are shown in Fig. [5ja). 
Due to the semi-metallic nature of graphene and the 
dense set of intra band transitions in the 0-5 eV energy 
region, a much denser k-point sampling is required to ob- 
tain a smooth absorption spectrum for this system. We 
used a 80 x 80 Monkhorst-Pack grid for both the ALDA 
and BSE calculations. 70 unoccupied bands were taken 
into account for the calculation of the response function, 
while 2 valence and 2 conduction band were included in 
the BSE Hamiltonian. The energy range below 1 eV is 
not shown in the figure since the excitations close to the 
Dirac point requires even denser k-points sampling. 

The LDA-ALDA spectrum (dashed-dotted line) shows 
absorption peaks at 3.9 and 5.6 eV originating from tran- 
sitions within the graphene and BN layer, respectively. 
It closely resembles a superposition of the spectra from 
freestanding graphene (not shown here) and BN sheets 
(dashed-dotted line in Fig. [3]) , with only a minor differ- 
ence of 0.1 eV in peak positions. Using GLLBSC-ALDA 
(dashed line), the two peaks shift up to 4.2 and 6.8 eV, 
respectively. The shift in the BN peak position is in ac- 
cordance with the shift in the BN gap in Fig. @] Note 
that the graphene peak energy of 4.2 eV is much lower 
than the 5.15 eV obtained from a previous Go Wo calcula- 
tion (without electron-hole interaction)^. We speculate 
that the deviation is due to an incorrect description of 
the slope of the graphene bands around the Dirac point 
where GLLBSC yields essentially the LDA result, see Fig. 
HJ Although the absolute absorption peak for graphene is 
underestimated, the excitonic effect is still well described 
using the BSE. With electron-hole pair interaction in- 
cluded (solid line), the graphene absorption peak at 4.2 
eV is redshifted by 0.6 eV, the same amount as was found 
in Ref. HU. The shift in the BN peak is, however, more 
striking. Upon adsorption of graphene, the BN exciton 



X (rir 2 ;r 3 r 4 ,cj) = P (r x r 2 ; r 3 r 4 ,u;) + / P (rir 2 ; 



peak shifts from 6.1 eV in Fig. [3]to 5.4 eV in Fig. [5] The 
reduction of the exciton energy of 0.7 eV is much smaller 
than the 1.98 eV reduction of the fundamental gap. This 
means that the exciton binding energy has been reduced 
from 1.9 eV in freestanding BN to 0.6 eV when adsorbed 
on graphene. Again, this is explained by the enhanced 
screening of the electron-hole pair provided by the elec- 
trons in graphene. The substrate induced screening of ex- 
citon binding energies was recently observed in GW-BSE 
calculations for molecules adsorbed on a metal surface* 6 - 



VI. CONCLUSIONS 

We have presented an implementation of the Bethe- 
Salpeter equation (BSE) which allows for the calculation 
of optical properties of materials with proper account 
of electron-hole interactions. Rather than following the 
standard approach where quasiparticle energies are ob- 
tained from the computationally costly GW method, we 
showed that excellent agreement with experimental ab- 
sorption spectra of a representative set of semiconductors 
and insulators, can be obtained by using single-particle 
energies from the GLLBSC functional. The latter yields 
very good fundamental gaps due to its explicit inclusion 
of the derivative discontinuity, and its computational cost 
is comparable to LDA. For a single layer of boron-nitride 
the fundamental gap and optical spectrum obtained with 
GLLBSC-BSE is very close to that of previous GW-BSE 
calculations. We showed that when BN is adsorbed on 
graphene, the fundamental gap is reduced by 2 eV. This 
reduction can be explained by image charge screening, 
and shows up in the GLLBSC calculation as a vanishing 
contribution from the derivative discontinuity. Finally, 
we found that the absoption spectrum of graphenc/BN 
interface is not simply a sum of the absorption spectra of 
the isolated layers, because the transition energies in BN 
become redshifted by up to almost 1 eV due to screening 
by the graphene electrons. 



Appendix A: Effective two-particle Hamiltonian 

To obtain an effective two-particle Hamiltonian de- 
scribing the optical excitations of the interacting electron 
system, we begin by considering the Bethe-Salpeter equa- 
tion (BSE) for the (retarded) four-point response func- 
tion, x 4P - Assuming a static electron- hole interaction 
kernel x 4P can be written 



,v)K (r 5 r 6 ;r 7 r 8 )x (r 7 r 8 ; r 3 r 4 , uj)dr 5 dr 6 dr 7 dr$ (Al) 
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In writing the above BSE equation we have made the sim- 
plifying, and for practical purposes essential, assumption 
that the electron-hole interaction kernel, K, is frequency 
independent. The quantity x 4P is an uncontracted ver- 
sion of the density response function, i.e. x( r ! r '; CJ ) = 
X 4P (rr; rV, uS) while P 4P is the four-point response func- 
tion for independent (but self-energy dressed) quasipar- 
ticles (QP). The kernel is given by K 4P = V — \W where 



V(rir 2 ;r 3 r 4 ) = 



1 



n - r 3 



-Mv x - r 2 )<5(r 3 - r 4 ) (A2) 



is the electron-hole exchange and 



T/F(rir 2 ;r 3 r 4 ) 



e ^r^r ,0) , 

dr d(ri - r 3 )S(r 2 - r 4 ) 



(A3) 



1-2 | 



is the statically screened direct electron-hole interaction. 

Assuming that the QP energies and wave functions can 
be described by an effective non-interacting Hamiltonian, 
Hqp , we can write the independent response function as 



P 4P (rir 2 ;r 3 r4,a;) = - ^ (/™k - /mk+q) 

q knm 

Ck ( r l ) V'mk+q (r 2 ) V'nk (r 3 ) V> mk+q (r 4 ) 
X , QP QP i • 

W + Vk- e mk+q + %T 1 



(A4) 



where the wave functions form an orthonormal set and 
the occupation factors are 1 or for occupied and empty 
states, respectively. 

The full four-point response function can also be ex- 
panded in the orthonormal basis of single-particle tran- 
sitions, ^ s (ri,r 2 ) = V£k( r l) 1 /'mk+q(r2), 



X 4P (rir 2 ;r 3 r 4 ,w) = ^^xss'(q,w) 

q SS' 

X Ck( r l)V' m k+q(r 2 )V'™'k' (l - 3)C l 'k'+q( r 4) 



(A5) 



As a consequence of the periodicity of the crystal lattice, 
all 4-point functions are diagonal in q. Note that the in- 
dices n, m, n' , m! must run over all bands, both occupied 
and unoccupied, in order to ensure that the two-particle 
basis is complete (it will, however, turn out that it is 
sufficient to consider only e-h and h-e transitions). 

The non-interacting response function is diagonal in 
the two-particle basis, 



fs 



uj - es + W) 



-Sss> 



(A6) 



where the occupation and transition energy for an elec- 
tron hole pair S is defined as 



(A7) 
(A8) 



fs — /nk /mk+q 
E-S = ^-nk ^mk+q 



The four point Bcthc-Salpeter equation Eq. (|Al|) in 
the two-particle basis corresponding to momentum trans- 



fer q becomes 

X 4 / 5 ,(q, W ) = P 4P .( q ,w) (A9) 

+ £ P 4P (q, u)Kf s „ (q, w) X 4P s , (q, W ) 
S" 

Expressions for the kernel matrix elements are given in 
Eqs. (PH and (fTT]). 

Substituting Eq. (|A6|) into Eq. (|A9|) and rearranging 
yields 



(A10) 



where the effective two-particle Hamiltonian, %, is de- 
fined as 



Hss>(q.,w) = e s 5 S s> + f s Kfg,(q,uj), 



(All) 



and I is an identity matrix with the same dimension as 
H. 

By dividing the matrices into 4x4 blocks correspond- 
ing to two-particle basis functions containing e-h, h-e, 
e-e, and h-h transitions, it follows that xWs' ^ s non-zero 
only within the 2x2 upper left block. For this reason 
we can reduce the problem by limiting the two-particle 
basis functions, ips, to the e-h and h-e states. Using the 
eigenstates and energies of the BSE Hamiltonian, 



H(q)A x (q) = E x (q)A x (q) 



(A12) 



we can construct the spectral representation of the resol- 
vent of the BSE Hamiltonian, 



^ Af(q)[Af,(q)]*jV- 4 (q) 

AA' 



u - E\(q) + ir) 



(A13) 

where N\\> (q) is the overlap matrix defined as 

iWq) = $>f(q)]Mf,(q) (AH) 



The BSE Hamiltonian (| A12I) is in general non-Hermitian 
as a matrix in the e-h and h-e basis. However, within 
the standard Tamm-Dancoff approximation, in which 
only the e-h transitions are considered (i.e. transitions 
with positive energies), H(q) becomes Hcrmitian and 
N\\>(q) = <5 A A'- 

Since the two-point response function, XGcfaiW), is 
obtained by Fourier transforming x 4P (rr; r'r', w), we con- 
clude from Eq. (|A5j) that 



Xgg' (q, V) = - J2 Xfs> ^s(G)n* s , (G') (A15) 

SS' 

where the charge density matrix, ns(G), is defined in Eq. 
©• 

Finally, the relation to the macroscopic dielectric func- 
tion Eq. ([TBJ is established using Eqs. (|A10|) and (|A13J) . 
together with the relation 



e(u) = 1 



47T 

w 



xoo(q^ o,w) 



(A16) 
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between the dielectric function and the irreducible re- 
sponse function, x.- As discussed in Sec. Ill CI the latter is 
obtained in place of x when the long range G = term 
excluded from the e-h exchange kernel in Eq. (|TU)) . 
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